home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1993 July / InfoMagic USENET CD-ROM July 1993.ISO / sources / misc / volume15 / ggems / part02 < prev    next >
Encoding:
Text File  |  1990-10-14  |  55.6 KB  |  2,120 lines

  1. Newsgroups: comp.sources.misc
  2. X-UNIX-From: craig@weedeater.math.yale.edu
  3. subject: v15i024: Graphics Gems, Part 2/5
  4. from: Craig Kolb <craig@weedeater.math.yale.edu>
  5. Sender: allbery@uunet.UU.NET (Brandon S. Allbery - comp.sources.misc)
  6.  
  7. Posting-number: Volume 15, Issue 24
  8. Submitted-by: Craig Kolb <craig@weedeater.math.yale.edu>
  9. Archive-name: ggems/part02
  10.  
  11. #! /bin/sh
  12. # This is a shell archive.  Remove anything before this line, then unpack
  13. # it by saving it into a file and typing "sh file".  To overwrite existing
  14. # files, type "sh file -c".  You can also feed this as standard input via
  15. # unshar, or by typing "sh <file", e.g..  If this archive is complete, you
  16. # will see the following message at the end:
  17. #        "End of archive 2 (of 5)."
  18. # Contents:  2DClip/bio.c 2DClip/clip.c BoxSphere.c Label.c Makefile
  19. #   MatrixPost.c Median.c OrderDither.c PntOnLine.c
  20. #   PolyScan/fancytest.c PolyScan/poly.c PolyScan/poly.h
  21. #   PolyScan/scantest.c Quaternions.c SeedFill.c SquareRoot.c
  22. #   Sturm/main.c TriPoints.c
  23. # Wrapped by craig@weedeater on Fri Oct 12 15:53:12 1990
  24. PATH=/bin:/usr/bin:/usr/ucb ; export PATH
  25. if test -f '2DClip/bio.c' -a "${1}" != "-c" ; then 
  26.   echo shar: Will not clobber existing file \"'2DClip/bio.c'\"
  27. else
  28. echo shar: Extracting \"'2DClip/bio.c'\" \(3146 characters\)
  29. sed "s/^X//" >'2DClip/bio.c' <<'END_OF_FILE'
  30. X
  31. X/*
  32. X * file  bio.c
  33. X *    contains the basic operations
  34. X *
  35. X */
  36. X#include    <stdio.h>
  37. X#include    "GraphicsGems.h"
  38. X#include    "line.h"
  39. X
  40. X/*
  41. X * def_contour
  42. X *
  43. X *    Purpose:
  44. X *    add a contour to the list
  45. X *    NOTE: coordinates must already be converted into longs!
  46. X *    
  47. X *    x    x coordinates of the end points of segments
  48. X *    y    y coordinates of the end points of segments
  49. X *    n    number of coordinate pairs
  50. X *    no    contour number (id), does no have to be unique!
  51. X *    type    type of clip operation desired CLIP_NORMAL means
  52. X *        clip everything inside the contour
  53. X */
  54. Xdef_contour(x, y, n, no, type)
  55. Xlong    x[], y[];
  56. Xint    n, no, type;
  57. X{
  58. X    short    i;
  59. X    long    dx1, dx2, dy1, dy2;
  60. X    long    minx, miny, maxx, maxy;
  61. X    CONTOUR    *cp;
  62. X    SEGMENT    *sp, *spp;
  63. X
  64. X    if((cp = CL) == (CONTOUR *)NULL) {
  65. X        cp = CL = NEWTYPE(CONTOUR);
  66. X    }
  67. X    else {
  68. X        while(cp->_next != (CONTOUR *)NULL)
  69. X            cp = cp->_next;
  70. X        i = cp->_no;
  71. X        cp = cp->_next = NEWTYPE(CONTOUR);
  72. X    }
  73. X
  74. X    cp->_next = (CONTOUR *)NULL;
  75. X    cp->_no = no;
  76. X    SET_ON(cp);
  77. X    if(type == CLIP_NORMAL)
  78. X        SET_INVERSE(cp);
  79. X    else
  80. X        SET_NORMAL(cp);
  81. X    minx = miny = 1000000;
  82. X    maxx = maxy = -1;
  83. X    for(i=0; i<n; i++) {
  84. X        if(i == 0) {
  85. X            cp->_s = sp = NEWTYPE(SEGMENT);
  86. X            sp->_from._x = x[0];
  87. X            sp->_from._y = y[0];
  88. X            sp->_to._x   = x[1];
  89. X            sp->_to._y   = y[1];
  90. X        }
  91. X        else {
  92. X        /*
  93. X         * if necessary stretch the contour
  94. X         * and skip the point
  95. X         */
  96. X            dx1 = sp->_to._x - sp->_from._x;
  97. X            dy1 = sp->_to._y - sp->_from._y;
  98. X            dx2 = x[(i == n-1) ? 0 : i+1] - sp->_to._x;
  99. X            dy2 = y[(i == n-1) ? 0 : i+1] - sp->_to._y;
  100. X            if(dy2*dx1 == dy1*dx2) {
  101. X                sp->_to._x = x[(i == n-1) ? 0 : i+1];
  102. X                sp->_to._y = y[(i == n-1) ? 0 : i+1];
  103. X            }
  104. X            else {
  105. X                spp = sp;
  106. X                sp = sp->_next =  NEWTYPE(SEGMENT);
  107. X                sp->_prev = spp;
  108. X                sp->_from._x = x[i];
  109. X                sp->_from._y = y[i];
  110. X                sp->_to._x = x[(i == n-1) ? 0 : i+1];
  111. X                sp->_to._y = y[(i == n-1) ? 0 : i+1];
  112. X            }
  113. X        }
  114. X
  115. X/*
  116. X * calculate the enclosing box
  117. X */
  118. X        if(x[i] < minx) 
  119. X            minx = x[i];
  120. X        if(x[i] > maxx)
  121. X            maxx = x[i];
  122. X        if(y[i] < miny)
  123. X            miny = y[i];
  124. X        if(y[i] > maxy)
  125. X            maxy = y[i];
  126. X            
  127. X    }
  128. X    cp->_minx = minx;
  129. X    cp->_maxx = maxx;
  130. X    cp->_miny = miny;
  131. X    cp->_maxy = maxy;
  132. X    sp->_next = cp->_s;
  133. X    cp->_s->_prev = sp;
  134. X    cp->_next = (CONTOUR *)NULL;
  135. X}
  136. X
  137. X/*
  138. X * get_contour_ptr
  139. X *
  140. X *    PURPOSE
  141. X *    get the pointer to a contour given its id
  142. X *    with multiple id's first fit algorithm is
  143. X *    used. Returns NULL in case of error.
  144. X *
  145. X *    no    id of contour
  146. X */
  147. XCONTOUR    *get_contour_ptr(no)
  148. Xint    no;
  149. X{
  150. X    CONTOUR    *cp;
  151. X
  152. X    if((cp = CL) == (CONTOUR *)NULL) 
  153. X        return((CONTOUR *)NULL);
  154. X    else {
  155. X        while(cp != (CONTOUR *)NULL) {
  156. X            if(cp->_no == no)
  157. X                return(cp);
  158. X            else
  159. X                cp = cp->_next;
  160. X        }
  161. X        return((CONTOUR *)NULL);
  162. X    }
  163. X}
  164. X
  165. X
  166. X/*
  167. X * del_contour
  168. X *
  169. X *    PURPOSE
  170. X *    delete contour(s) from the list with id
  171. X *    no
  172. X */
  173. Xdel_contour(no)
  174. Xint    no;
  175. X{
  176. X    CONTOUR    *cp, *cpp;
  177. X    CONTOUR    *qp = (CONTOUR *)NULL;
  178. X    SEGMENT    *sp, *spp;
  179. X
  180. X    if((cp = CL) == (CONTOUR *)NULL)
  181. X        return;
  182. X    while(cp != (CONTOUR *)NULL) {
  183. X        if(cp->_no == no) {
  184. X            sp = cp->_s;
  185. X            do {
  186. X                spp = sp->_next;
  187. X                free(sp);
  188. X                sp = spp;
  189. X            } while(sp != cp->_s);
  190. X            cpp = cp->_next;
  191. X            free(cp);
  192. X            if(qp)
  193. X                qp->_next = cpp;
  194. X            else
  195. X                CL = cpp;
  196. X            cp = cpp;
  197. X        } 
  198. X        else {
  199. X            qp = cp;
  200. X            cp = cp->_next;
  201. X        }
  202. X    }
  203. X}
  204. X
  205. END_OF_FILE
  206. if test 3146 -ne `wc -c <'2DClip/bio.c'`; then
  207.     echo shar: \"'2DClip/bio.c'\" unpacked with wrong size!
  208. fi
  209. # end of '2DClip/bio.c'
  210. fi
  211. if test -f '2DClip/clip.c' -a "${1}" != "-c" ; then 
  212.   echo shar: Will not clobber existing file \"'2DClip/clip.c'\"
  213. else
  214. echo shar: Extracting \"'2DClip/clip.c'\" \(2904 characters\)
  215. sed "s/^X//" >'2DClip/clip.c' <<'END_OF_FILE'
  216. X/*
  217. X * file clip.c
  218. X *    contains the actual clipping routines
  219. X */
  220. X#include    <stdio.h>
  221. X#include    "GraphicsGems.h"
  222. X#include    "line.h"
  223. X
  224. X/*
  225. X * vis_vector
  226. X *
  227. X *    PURPOSE
  228. X *    actual user interface. Draws a clipped line
  229. X *    NOTE: coordinates are given in converted LONGS!
  230. X *
  231. X *    xf, yf    from coordinates of vector to be drawn
  232. X *    xt, yt    to coordinates of vector to be drawn
  233. X */
  234. Xvis_vector(xf, yf, xt, yt)
  235. Xlong    xf, yf, xt, yt;
  236. X{
  237. X    SEGMENT    l;
  238. X
  239. X    if(xf == xt && yf == yt)
  240. X        return;
  241. X    l._from._x = xf;
  242. X    l._from._y = yf;
  243. X    l._to._x   = xt;
  244. X    l._to._y   = yt;
  245. X/*
  246. X * start at top of list
  247. X */
  248. X    clip(CL, &l);
  249. X}
  250. X
  251. X/*
  252. X * clip
  253. X *
  254. X *    PURPOSE
  255. X *
  256. X *    p    pointer to polygon
  257. X *    l    pointer to line segment
  258. X */
  259. Xclip(p, l)
  260. XCONTOUR    *p;
  261. XSEGMENT    *l;
  262. X{
  263. X    SEGMENT    *sp, ss;
  264. X    CLIST    *sol;
  265. X    POINT    pt;
  266. X    boolean    up, delay, inside, p_inside(), disjunct();
  267. X    int    i;
  268. X    short    nsol, nsmax = 2;
  269. X
  270. X
  271. X/*
  272. X * list exhausted do what you like
  273. X * we want to plot
  274. X */
  275. X    if(p == (CONTOUR *)NULL) {
  276. X        move((l->_from._x), (l->_from._y));
  277. X        cont((l->_to._x), (l->_to._y));
  278. X        return;
  279. X    }
  280. X/*
  281. X * polygon is switched off
  282. X * take next one
  283. X */
  284. X    if(!IS_ON(p)) {
  285. X        clip(p->_next, l);
  286. X        return;
  287. X    }
  288. X/*
  289. X * comparison on basis of the
  290. X * enclosing rectangle
  291. X */
  292. X    if(disjunct(p, l)) {
  293. X        if(!IS_NORMAL(p)) {
  294. X            clip(p->_next, l);
  295. X        }
  296. X        return;
  297. X    }
  298. X/*
  299. X * calculate possible intersections
  300. X */
  301. X    sol = (CLIST *) calloc(2, sizeof(CLIST));
  302. X    sol[0]._p._x = l->_from._x;
  303. X    sol[0]._p._y = l->_from._y;
  304. X    sol[0]._type = STD;
  305. X    sol[1]._p._x = l->_to._x;
  306. X    sol[1]._p._y = l->_to._y;
  307. X    sol[1]._type = STD;
  308. X    nsol = 2;
  309. X    cross_calc(p, l, &sol, &nsol, nsmax);
  310. X    pt._x = sol[0]._p._x;
  311. X    pt._y = sol[0]._p._y;
  312. X
  313. X/*
  314. X * determine status of first point
  315. X */
  316. X    inside = p_inside(p, &pt);
  317. X    if((!inside && IS_NORMAL(p)) || (inside && !IS_NORMAL(p)))
  318. X        up = TRUE; 
  319. X    else
  320. X        up = FALSE;
  321. X    delay = FALSE;
  322. X/*
  323. X * process list of intersections
  324. X */
  325. X    for(i=1; i<nsol; i++) {
  326. X        if(!up) {
  327. X            ss._from._x = sol[i-1]._p._x;
  328. X            ss._from._y = sol[i-1]._p._y;
  329. X            ss._to._x = sol[i]._p._x;
  330. X            ss._to._y = sol[i]._p._y;
  331. X            clip(p->_next, &ss);
  332. X        }
  333. X        if(!delay) {
  334. X            if(sol[i]._type != DELAY)
  335. X                up = (up) ? FALSE : TRUE;
  336. X            else
  337. X                delay = TRUE;
  338. X        }
  339. X        else {
  340. X            up = (up) ? FALSE : TRUE;
  341. X            delay = FALSE;
  342. X        }
  343. X    }
  344. X    free(sol);
  345. X}
  346. X
  347. X/*
  348. X * disjunct
  349. X *
  350. X *    PURPOSE
  351. X *    determine if the box enclosing the polygon 
  352. X *    stored in p and the box enclosing the line 
  353. X *    segment stored in l are disjunct.
  354. X *    Return TRUE if disjunct else FALSE
  355. X *
  356. X *    p    points to the polygon structure
  357. X *    l    points to the linesegment structure    
  358. X *
  359. X */
  360. Xboolean    disjunct(p, l)
  361. XCONTOUR    *p;
  362. XSEGMENT    *l;
  363. X
  364. X{
  365. X    if((max(l->_from._x, l->_to._x) < p->_minx) ||
  366. X       (min(l->_from._x, l->_to._x) > p->_maxx) ||
  367. X           (max(l->_from._y, l->_to._y) < p->_miny) ||
  368. X       (min(l->_from._y, l->_to._y) > p->_maxy)   )
  369. X        return(TRUE);
  370. X    else
  371. X        return(FALSE);
  372. X}
  373. X
  374. X#define DEBUG
  375. X#ifdef DEBUG
  376. Xmove(x, y)
  377. Xlong    x, y;
  378. X{
  379. X    printf("(%d,%d) ->", x, y);
  380. X}
  381. X
  382. Xcont(x, y)
  383. Xlong    x, y;
  384. X{
  385. X    printf("(%d,%d)\n", x, y);
  386. X}
  387. X
  388. X#endif
  389. X
  390. X
  391. X
  392. X
  393. END_OF_FILE
  394. if test 2904 -ne `wc -c <'2DClip/clip.c'`; then
  395.     echo shar: \"'2DClip/clip.c'\" unpacked with wrong size!
  396. fi
  397. # end of '2DClip/clip.c'
  398. fi
  399. if test -f 'BoxSphere.c' -a "${1}" != "-c" ; then 
  400.   echo shar: Will not clobber existing file \"'BoxSphere.c'\"
  401. else
  402. echo shar: Extracting \"'BoxSphere.c'\" \(3760 characters\)
  403. sed "s/^X//" >'BoxSphere.c' <<'END_OF_FILE'
  404. X/*
  405. XA Simple Method for Box-Sphere Intersection Testing
  406. Xby Jim Arvo
  407. Xfrom "Graphics Gems", Academic Press, 1990
  408. X*/
  409. X
  410. X#include "GraphicsGems.h"
  411. X
  412. X/*
  413. X *  This routine tests for intersection between an n-dimensional             
  414. X *  axis-aligned box and an n-dimensional sphere.  The mode argument       
  415. X *  indicates whether the objects are to be regarded as surfaces or        
  416. X *  solids.  The values are:                                               
  417. X *                                                                         
  418. X *     mode                                                                
  419. X *                                                                         
  420. X *       0   Hollow Box, Hollow Sphere                                     
  421. X *       1   Hollow Box, Solid  Sphere
  422. X *       2   Solid  Box, Hollow Sphere                                     
  423. X *       3   Solid  Box, Solid  Sphere                                     
  424. X *                                                                         
  425. X*/
  426. Xint Box_Sphere_Intersect( n, Bmin, Bmax, C, r, mode )
  427. Xint    n;       /* The dimension of the space.           */
  428. Xfloat  Bmin[];  /* The minium of the box for each axis.  */
  429. Xfloat  Bmax[];  /* The maximum of the box for each axis. */
  430. Xfloat  C[];     /* The sphere center in n-space.         */
  431. Xfloat  r;       /* The radius of the sphere.             */
  432. Xint    mode;    /* Selects hollow or solid.              */
  433. X    {
  434. X    float  a, b;
  435. X    float  dmin, dmax;
  436. X    float  r2 = SQR( r );
  437. X    int    i, face;
  438. X
  439. X
  440. X    switch( mode )
  441. X        {
  442. X        case 0: /* Hollow Box - Hollow Sphere */
  443. X            dmin = 0;
  444. X            dmax = 0;
  445. X            face = FALSE;
  446. X            for( i = 0; i < n; i++ ) {
  447. X                a = SQR( C[i] - Bmin[i] );
  448. X                b = SQR( C[i] - Bmax[i] );
  449. X                dmax += MAX( a, b );
  450. X                if( C[i] < Bmin[i] ) {
  451. X                    face = TRUE;
  452. X                    dmin += a;
  453. X                    }
  454. X                else if( C[i] > Bmax[i] ) {
  455. X                    face = TRUE;
  456. X                    dmin += b;
  457. X                    }
  458. X                else if( MIN( a, b ) <= r2 ) face = TRUE;
  459. X                }
  460. X            if(face && ( dmin <= r2 ) && ( r2 <= dmax)) return(TRUE);
  461. X            break;
  462. X
  463. X        case 1: /* Hollow Box - Solid Sphere */
  464. X            dmin = 0;
  465. X            face = FALSE;
  466. X            for( i = 0; i < n; i++ ) {
  467. X                if( C[i] < Bmin[i] ) {
  468. X                    face = TRUE;
  469. X                    dmin += SQR( C[i] - Bmin[i] );
  470. X                    }
  471. X                else if( C[i] > Bmax[i] ) {
  472. X                    face = TRUE;
  473. X                    dmin += SQR( C[i] - Bmax[i] );     
  474. X                    }
  475. X                else if( C[i] - Bmin[i] <= r ) face = TRUE;
  476. X                else if( Bmax[i] - C[i] <= r ) face = TRUE;
  477. X                }
  478. X            if( face && ( dmin <= r2 ) ) return( TRUE );
  479. X            break;
  480. X
  481. X
  482. X        case 2: /* Solid Box - Hollow Sphere */
  483. X            dmax = 0;
  484. X            dmin = 0;
  485. X            for( i = 0; i < n; i++ ) {
  486. X                a = SQR( C[i] - Bmin[i] );
  487. X                b = SQR( C[i] - Bmax[i] );
  488. X                dmax += MAX( a, b );
  489. X                if( C[i] < Bmin[i] ) dmin += a; else
  490. X                if( C[i] > Bmax[i] ) dmin += b;
  491. X                }
  492. X            if( dmin <= r2 && r2 <= dmax ) return( TRUE );
  493. X            break;
  494. X
  495. X        case 3: /* Solid Box - Solid Sphere */
  496. X            dmin = 0;
  497. X            for( i = 0; i < n; i++ ) {
  498. X                if( C[i] < Bmin[i] ) dmin += SQR(C[i] - Bmin[i] ); else
  499. X                if( C[i] > Bmax[i] ) dmin += SQR( C[i] - Bmax[i] );     
  500. X                }
  501. X            if( dmin <= r2 ) return( TRUE );
  502. X            break;
  503. X  
  504. X        } /* end switch */
  505. X
  506. X    return( FALSE );
  507. X    } 
  508. END_OF_FILE
  509. if test 3760 -ne `wc -c <'BoxSphere.c'`; then
  510.     echo shar: \"'BoxSphere.c'\" unpacked with wrong size!
  511. fi
  512. # end of 'BoxSphere.c'
  513. fi
  514. if test -f 'Label.c' -a "${1}" != "-c" ; then 
  515.   echo shar: Will not clobber existing file \"'Label.c'\"
  516. else
  517. echo shar: Extracting \"'Label.c'\" \(2291 characters\)
  518. sed "s/^X//" >'Label.c' <<'END_OF_FILE'
  519. X/*
  520. X * Nice Numbers for Graph Labels
  521. X * by Paul Heckbert
  522. X * from "Graphics Gems", Academic Press, 1990
  523. X */
  524. X
  525. X/*
  526. X * label.c: demonstrate nice graph labeling
  527. X *
  528. X * Paul Heckbert    2 Dec 88
  529. X */
  530. X
  531. X#include <stdio.h>
  532. X#include <math.h>
  533. X#include "GraphicsGems.h"
  534. X
  535. Xdouble nicenum();
  536. X
  537. X/* expt(a,n)=a^n for integer n */
  538. X
  539. X#ifdef POW_NOT_TRUSTWORTHY
  540. X/* if roundoff errors in pow cause problems, use this: */
  541. X
  542. Xdouble expt(a, n)
  543. Xdouble a;
  544. Xregister int n;
  545. X{
  546. X    double x;
  547. X
  548. X    x = 1.;
  549. X    if (n>0) for (; n>0; n--) x *= a;
  550. X    else for (; n<0; n++) x /= a;
  551. X    return x;
  552. X}
  553. X
  554. X#else
  555. X#   define expt(a, n) pow(a, (double)(n))
  556. X#endif
  557. X
  558. X#define NTICK 5            /* desired number of tick marks */
  559. X
  560. Xmain(ac, av)
  561. Xint ac;
  562. Xchar **av;
  563. X{
  564. X    double min, max;
  565. X
  566. X    if (ac!=3) {
  567. X    fprintf(stderr, "Usage: label <min> <max>\n");
  568. X    exit(1);
  569. X    }
  570. X    min = atof(av[1]);
  571. X    max = atof(av[2]);
  572. X
  573. X    loose_label(min, max);
  574. X}
  575. X
  576. X/*
  577. X * loose_label: demonstrate loose labeling of data range from min to max.
  578. X * (tight method is similar)
  579. X */
  580. X
  581. Xloose_label(min, max)
  582. Xdouble min, max;
  583. X{
  584. X    char str[6], temp[20];
  585. X    int nfrac;
  586. X    double d;                /* tick mark spacing */
  587. X    double graphmin, graphmax;        /* graph range min and max */
  588. X    double range, x;
  589. X
  590. X    /* we expect min!=max */
  591. X    range = nicenum(max-min, 0);
  592. X    d = nicenum(range/(NTICK-1), 1);
  593. X    graphmin = floor(min/d)*d;
  594. X    graphmax = ceil(max/d)*d;
  595. X    nfrac = MAX(-floor(log10(d)), 0);    /* # of fractional digits to show */
  596. X    sprintf(str, "%%.%df", nfrac);    /* simplest axis labels */
  597. X
  598. X    printf("graphmin=%g graphmax=%g increment=%g\n", graphmin, graphmax, d);
  599. X    for (x=graphmin; x<graphmax+.5*d; x+=d) {
  600. X    sprintf(temp, str, x);
  601. X    printf("(%s)\n", temp);
  602. X    }
  603. X}
  604. X
  605. X/*
  606. X * nicenum: find a "nice" number approximately equal to x.
  607. X * Round the number if round=1, take ceiling if round=0
  608. X */
  609. X
  610. Xstatic double nicenum(x, round)
  611. Xdouble x;
  612. Xint round;
  613. X{
  614. X    int exp;                /* exponent of x */
  615. X    double f;                /* fractional part of x */
  616. X    double nf;                /* nice, rounded fraction */
  617. X
  618. X    exp = floor(log10(x));
  619. X    f = x/expt(10., exp);        /* between 1 and 10 */
  620. X    if (round)
  621. X    if (f<1.5) nf = 1.;
  622. X    else if (f<3.) nf = 2.;
  623. X    else if (f<7.) nf = 5.;
  624. X    else nf = 10.;
  625. X    else
  626. X    if (f<=1.) nf = 1.;
  627. X    else if (f<=2.) nf = 2.;
  628. X    else if (f<=5.) nf = 5.;
  629. X    else nf = 10.;
  630. X    return nf*expt(10., exp);
  631. X}
  632. END_OF_FILE
  633. if test 2291 -ne `wc -c <'Label.c'`; then
  634.     echo shar: \"'Label.c'\" unpacked with wrong size!
  635. fi
  636. # end of 'Label.c'
  637. fi
  638. if test -f 'Makefile' -a "${1}" != "-c" ; then 
  639.   echo shar: Will not clobber existing file \"'Makefile'\"
  640. else
  641. echo shar: Extracting \"'Makefile'\" \(2894 characters\)
  642. sed "s/^X//" >'Makefile' <<'END_OF_FILE'
  643. X#
  644. X# Makefile for Graphics Gems source
  645. X#
  646. X# Craig Kolb, 8/90
  647. X#
  648. X# This make file will build "gemslib.a" and a number of executables.
  649. X# Gemslib is built solely for debugging purposes -- it is not intended
  650. X# to be used as a library.
  651. X#
  652. X# Note that some of the gems need additional macros, functions, tables
  653. X# driving routines, etc. before they will compile or run properly.
  654. X# These include:
  655. X#
  656. X# AALines
  657. X#    Needs variables set in the Makefile.
  658. X# Dissolve
  659. X#    Needs a table filled.
  660. X# FitCurves
  661. X#    Needs a functioning DrawBezierCurve().
  662. X# MatrixOrtho
  663. X#    Needs a number of matrix routines.
  664. X# PolyScan
  665. X#    Needs driving routines and other additional code.
  666. X# RayPolygon
  667. X#    Needs surrounding function structure, declaration of
  668. X#    variable types, etc.
  669. X
  670. X#
  671. X# C compiler flags
  672. X#
  673. XCFLAGS = -g
  674. X#
  675. X# Location of Graphics Gems library
  676. X#
  677. XLIBFILE = gemslib.a
  678. X
  679. X#
  680. X# Graphics Gems Vector Library
  681. X#
  682. XVECLIB = GGVecLib.o
  683. X
  684. XMFLAGS = "LIBFILE = ../$(LIBFILE)" "GENCFLAGS = $(CFLAGS)"
  685. X
  686. XFTPDIR = /usr/spool/uucppublic/ftp/pub/GraphicsGems/src
  687. XZOOFILE = Gems.zoo
  688. X
  689. XSHELL = /bin/sh
  690. X
  691. XOFILES =    PntOnLine.o ViewTrans.o AAPolyScan.o Albers.o \
  692. X        Interleave.o BoundSphere.o BoxSphere.o CircleRect.o \
  693. X        ConcaveScan.o Roots3And4.o Dissolve.o DigitalLine.o \
  694. X        FastJitter.o FixedTrig.o HSLtoRGB.o HypotApprox.o LineEdge.o \
  695. X        MatrixInvert.o MatrixPost.o Median.o PixelInteger.o \
  696. X        TriPoints.o Quaternions.o RGBTo4Bits.o RayBox.o \
  697. X        SeedFill.o SquareRoot.o DoubleLine.o TransBox.o
  698. X
  699. XDIRS = 2DClip PolyScan Sturm
  700. X
  701. XALL =    Hash3D FitCurves Forms NearestPoint Label \
  702. X    OrderDither BinRec $(LIBFILE)
  703. X
  704. Xall: $(ALL)
  705. X    @for d in $(DIRS) ; do \
  706. X        (cd $$d ; $(MAKE) $(MFLAGS)) ;\
  707. X    done
  708. X
  709. X$(LIBFILE): $(OFILES) $(VECLIB)
  710. X    ar rcs $(LIBFILE) $(OFILES) $(VECLIB)
  711. X
  712. XHash3D: Hash3D.o
  713. X    $(CC) $(CFLAGS) -o $@ Hash3D.o
  714. X
  715. XFitCurves: FitCurves.o $(VECLIB)
  716. X    $(CC) $(CFLAGS) -o $@ FitCurves.o $(VECLIB) -lm
  717. X
  718. XForms: Forms.o
  719. X    $(CC) $(CFLAGS) -o $@ Forms.o
  720. X
  721. XNearestPoint: NearestPoint.o $(VECLIB)
  722. X    $(CC) $(CFLAGS) -o $@ NearestPoint.o $(VECLIB) -lm
  723. X
  724. XLabel: Label.o
  725. X    $(CC) $(CFLAGS) -o $@ Label.o -lm
  726. X
  727. XOrderDither: OrderDither.o
  728. X    $(CC) $(CFLAGS) -o $@ OrderDither.o
  729. X
  730. XBinRec: BinRec.o
  731. X    $(CC) $(CFLAGS) -o $@ BinRec.o
  732. X
  733. Xftpreadme:
  734. X    /bin/rm -f $(FTPDIR)/README
  735. X    sed -e "s/__DATE__/`date`/g" < README.ftp > $(FTPDIR)/README
  736. X
  737. Xreadme:
  738. X    /bin/rm -rf README
  739. X    sed -e "s/__DATE__/`date`/g" < README.dist > README
  740. X    
  741. X
  742. Xclean:
  743. X    @for d in $(DIRS) ; do \
  744. X        (cd $$d ; $(MAKE) $(MFLAGS) clean) ;\
  745. X    done
  746. X    /bin/rm -f $(OFILES) $(VECLIB)
  747. X    /bin/rm -f Hash3D.o FitCurves.o \
  748. X        Forms.o NearestPoint.o Label.o OrderDither.o BinRec.o \
  749. X        Hash3D FitCurves Forms NearestPoint Label OrderDither BinRec \
  750. X        bugs a.out core Part??.Z Part?? $(ZOOFILE)
  751. X
  752. Xkit: readme
  753. X    /bin/rm -f Part??.Z Part??
  754. X    (makekit -iPACKING_LIST -oMANIFEST)
  755. X
  756. Xzoo: readme
  757. X    /bin/rm -f Gems.zoo
  758. X    cut -d" " -f2 PACKING_LIST | zoo aI $(ZOOFILE)
  759. X
  760. Xftp: kit zoo ftpreadme
  761. X    compress Part??
  762. X    /bin/cp Part??.Z $(ZOOFILE) $(FTPDIR)
  763. X
  764. X$(ALL): GraphicsGems.h
  765. END_OF_FILE
  766. if test 2894 -ne `wc -c <'Makefile'`; then
  767.     echo shar: \"'Makefile'\" unpacked with wrong size!
  768. fi
  769. # end of 'Makefile'
  770. fi
  771. if test -f 'MatrixPost.c' -a "${1}" != "-c" ; then 
  772.   echo shar: Will not clobber existing file \"'MatrixPost.c'\"
  773. else
  774. echo shar: Extracting \"'MatrixPost.c'\" \(3390 characters\)
  775. sed "s/^X//" >'MatrixPost.c' <<'END_OF_FILE'
  776. X/*
  777. XEfficient Post-Concatenation of Transformation Matrics
  778. Xby Joseph M. Cychosz
  779. Xfrom "Graphics Gems", Academic Press, 1990
  780. X*/
  781. X
  782. X#include    <math.h>
  783. X#include    "GraphicsGems.h"
  784. X
  785. X
  786. X/*       M4xform.c - Basic 4x4 transformation package.
  787. X *                                     
  788. X *    Description:                             
  789. X *        M4xform.c contains a collection of routines used to perform  
  790. X *        direct post-concatenated transformation operations.
  791. X *         
  792. X *                                     
  793. X *    Contents:                             
  794. X *        M4RotateX        Post-concatenate a x-axis rotation.     
  795. X *        M4RotateY        Post-concatenate a y-axis rotation.     
  796. X *        M4RotateZ        Post-concatenate a z-axis rotation.     
  797. X *        M4Scale            Post-concatenate a scaling.         
  798. X *        M4Translate        Post-concatenate a translation.         
  799. X *        M4ZPerspective    Post-concatenate a z-axis perspective     
  800. X *                        transformation.                 
  801. X *                                     
  802. X *    Externals:                             
  803. X *        cos, sin.                                             
  804. X */ 
  805. X     
  806. X
  807. X
  808. X/*      M4RotateX - Post-concatenate a x-axis rotation matrix.  */
  809. X
  810. XMatrix4    *M4RotateX    (m,a)
  811. X        Matrix4    *m;            /* Current transformation matrix*/
  812. X        double    a;            /* Rotation angle        */
  813. X{
  814. X        double    c, s;
  815. X        double    t;
  816. X        int        i;
  817. X
  818. X
  819. X        c = cos (a);    s = sin (a);
  820. X
  821. X        for (i = 0 ; i < 4 ; i++) {
  822. X             t = m->element[i][1];
  823. X             m->element[i][1] = t*c - m->element[i][2]*s;
  824. X             m->element[i][2] = t*s + m->element[i][2]*c;
  825. X        }
  826. X        return (m);
  827. X}
  828. X
  829. X
  830. X/*     M4RotateY - Post-concatenate a y-axis rotation matrix.  */
  831. X
  832. XMatrix4        *M4RotateY    (m,a)
  833. X        Matrix4    *m;            /* Current transformation matrix*/
  834. X        double    a;            /* Rotation angle        */
  835. X{
  836. X        double    c, s;
  837. X        double    t;
  838. X        int        i;
  839. X
  840. X        c = cos (a);    s = sin (a);
  841. X
  842. X        for (i = 0 ; i < 4 ; i++) {
  843. X             t = m->element[i][0];
  844. X             m->element[i][0] = t*c + m->element[i][2]*s;
  845. X             m->element[i][2] = m->element[i][2]*c - t*s;
  846. X        }
  847. X        return (m);
  848. X}
  849. X
  850. X
  851. X/*       M4RotateZ - Post-concatenate a z-axis rotation matrix.   */
  852. X
  853. XMatrix4    *M4RotateZ    (m,a)
  854. X        Matrix4    *m;            /* Current transformation matrix*/
  855. X        double    a;            /* Rotation angle        */
  856. X{
  857. X        double    c, s;
  858. X        double    t;
  859. X        int        i;
  860. X
  861. X        c = cos (a);    s = sin (a);
  862. X
  863. X        for (i = 0 ; i < 4 ; i++) {
  864. X             t = m->element[i][0];
  865. X             m->element[i][0] = t*c - m->element[i][1]*s;
  866. X             m->element[i][1] = t*s + m->element[i][1]*c;
  867. X        }
  868. X        return (m);
  869. X}
  870. X
  871. X
  872. X
  873. X/*       M4Scale    - Post-concatenate a scaling.   */
  874. X
  875. XMatrix4    *M4Scale    (m,sx,sy,sz)
  876. X        Matrix4    *m;            /* Current transformation matrix */
  877. X        double    sx, sy, sz;    /* Scale factors about x,y,z */
  878. X{
  879. X        int        i;
  880. X
  881. X        for (i = 0 ; i < 4 ; i++) {
  882. X             m->element[i][0] *= sx;
  883. X             m->element[i][1] *= sy;
  884. X             m->element[i][2] *= sz;
  885. X        }
  886. X        return (m);
  887. X}
  888. X
  889. X
  890. X/*       M4Translate - Post-concatenate a translation.   */
  891. X
  892. XMatrix4    *M4Translate    (m,tx,ty,tz)
  893. X        Matrix4    *m;            /* Current transformation matrix */
  894. X        double    tx, ty, tz;    /* Translation distance */
  895. X{
  896. X        int        i;
  897. X
  898. X        for (i = 0 ; i < 4 ; i++) {
  899. X             m->element[i][0] += m->element[i][3]*tx;
  900. X             m->element[i][1] += m->element[i][3]*ty;
  901. X             m->element[i][2] += m->element[i][3]*tz;
  902. X        }
  903. X        return (m);
  904. X}
  905. X
  906. X
  907. X        /* M4ZPerspective  Post-concatenate a perspective       */        /*transformation.                                    */
  908. X        /*                                                    */
  909. X        /* Perspective is along the z-axis with the eye at +z.   */
  910. X
  911. XMatrix4    *M4ZPerspective    (m,d)
  912. X        Matrix4    *m;            /* Current transformation matrix  */
  913. X        double    d;            /* Perspective distance        */
  914. X{
  915. X        int        i;
  916. X        double    f = 1. / d;
  917. X
  918. X        for (i = 0 ; i < 4 ; i++) {
  919. X             m->element[i][3] += m->element[i][2]*f;
  920. X             m->element[i][2]  = 0.;
  921. X        }
  922. X        return (m);
  923. X}
  924. X    
  925. X
  926. X
  927. END_OF_FILE
  928. if test 3390 -ne `wc -c <'MatrixPost.c'`; then
  929.     echo shar: \"'MatrixPost.c'\" unpacked with wrong size!
  930. fi
  931. # end of 'MatrixPost.c'
  932. fi
  933. if test -f 'Median.c' -a "${1}" != "-c" ; then 
  934.   echo shar: Will not clobber existing file \"'Median.c'\"
  935. else
  936. echo shar: Extracting \"'Median.c'\" \(2941 characters\)
  937. sed "s/^X//" >'Median.c' <<'END_OF_FILE'
  938. X/*
  939. XMedian Finding on a 3-by-3 Grid
  940. Xby Alan Paeth
  941. Xfrom "Graphics Gems", Academic Press, 1990
  942. X*/
  943. X
  944. X#define s2(a,b) {register int t; if ((t=b-a)<0) {a+=t; b-=t;}}
  945. X#define mn3(a,b,c) s2(a,b); s2(a,c);
  946. X#define mx3(a,b,c) s2(b,c); s2(a,c);
  947. X#define mnmx3(a,b,c) mx3(a,b,c); s2(a,b);
  948. X#define mnmx4(a,b,c,d) s2(a,b); s2(c,d); s2(a,c); s2(b,d);
  949. X#define mnmx5(a,b,c,d,e) s2(a,b); s2(c,d); mn3(a,c,e); mx3(b,d,e);
  950. X#define mnmx6(a,b,c,d,e,f) s2(a,d); s2(b,e); s2(c,f);\
  951. X                            mn3(a,b,c); mx3(d,e,f);
  952. Xmed3x3(b1, b2, b3)
  953. X    int *b1, *b2, *b3;
  954. X/*
  955. X * Find median on a 3x3 input box of integers.
  956. X * b1, b2, b3 are pointers to the left-hand edge of three
  957. X * parallel scan-lines to form a 3x3 spatial median.
  958. X * Rewriting b2 and b3 as b1 yields code which forms median
  959. X * on input presented as a linear array of nine elements.
  960. X */
  961. X    {
  962. X    register int r1, r2, r3, r4, r5, r6;
  963. X    r1 = *b1++; r2 = *b1++; r3 = *b1++;
  964. X    r4 = *b2++; r5 = *b2++; r6 = *b2++;
  965. X    mnmx6(r1, r2, r3, r4, r5, r6);
  966. X    r1 = *b3++;
  967. X    mnmx5(r1, r2, r3, r4, r5);
  968. X    r1 = *b3++;
  969. X    mnmx4(r1, r2, r3, r4);
  970. X    r1 = *b3++;
  971. X    mnmx3(r1, r2, r3);
  972. X    return(r2);
  973. X    }
  974. X
  975. X
  976. X/* t2(i,j) transposes elements in A[] such that A[i] <= A[j] */
  977. X
  978. X#define t2(i, j) s2(A[i-1], A[j-1])
  979. X
  980. X
  981. Xint median25(A)
  982. X    int A[25];
  983. X    {
  984. X/*
  985. X * median25(A) partitions the array A[0..24] such that element
  986. X * A[12] is the median and subarrays A[0..11] and A[13..24]
  987. X * are partitions containing elements of smaller and larger
  988. X * value (rank), respectively.
  989. X *
  990. X * The exchange table lists element indices on the range 1..25,
  991. X * this accounts for the "-1" offsets in the macro t2 and in
  992. X * the final return value used to adjust subscripts to C-code
  993. X * conventions (array indices begin at zero).
  994. X */
  995. X    t2( 1, 2); t2( 4, 5); t2( 3, 5); t2( 3, 4); t2( 7, 8);
  996. X    t2( 6, 8); t2( 6, 7); t2(10,11); t2( 9,11); t2( 9,10);
  997. X    t2(13,14); t2(12,14); t2(12,13); t2(16,17); t2(15,17);
  998. X    t2(15,16); t2(19,20); t2(18,20); t2(18,19); t2(22,23);
  999. X    t2(21,23); t2(21,22); t2(24,25); t2( 3, 6); t2( 4, 7);
  1000. X    t2( 1, 7); t2( 1, 4); t2( 5, 8); t2( 2, 8); t2( 2, 5);
  1001. X    t2(12,15); t2( 9,15); t2( 9,12); t2(13,16); t2(10,16);
  1002. X    t2(10,13); t2(14,17); t2(11,17); t2(11,14); t2(21,24);
  1003. X    t2(18,24); t2(18,21); t2(22,25); t2(19,25); t2(19,22);
  1004. X    t2(20,23); t2( 9,18); t2(10,19); t2( 1,19); t2( 1,10);
  1005. X    t2(11,20); t2( 2,20); t2( 2,11); t2(12,21); t2( 3,21);
  1006. X    t2( 3,12); t2(13,22); t2( 4,22); t2( 4,13); t2(14,23);
  1007. X    t2( 5,23); t2( 5,14); t2(15,24); t2( 6,24); t2( 6,15);
  1008. X    t2(16,25); t2( 7,25); t2( 7,16); t2( 8,17); t2( 8,20);
  1009. X    t2(14,22); t2(16,24); t2( 8,14); t2( 8,16); t2( 2,10);
  1010. X    t2( 4,12); t2( 6,18); t2(12,18); t2(10,18); t2( 5,11);
  1011. X    t2( 7,13); t2( 8,15); t2( 5, 7); t2( 5, 8); t2(13,15);
  1012. X    t2(11,15); t2( 7, 8); t2(11,13); t2( 7,11); t2( 7,18);
  1013. X    t2(13,18); t2( 8,18); t2( 8,11); t2(13,19); t2( 8,13);
  1014. X    t2(11,19); t2(13,21); t2(11,21); t2(11,13);
  1015. X    return(A[13-1]);
  1016. X    }
  1017. END_OF_FILE
  1018. if test 2941 -ne `wc -c <'Median.c'`; then
  1019.     echo shar: \"'Median.c'\" unpacked with wrong size!
  1020. fi
  1021. # end of 'Median.c'
  1022. fi
  1023. if test -f 'OrderDither.c' -a "${1}" != "-c" ; then 
  1024.   echo shar: Will not clobber existing file \"'OrderDither.c'\"
  1025. else
  1026. echo shar: Extracting \"'OrderDither.c'\" \(2483 characters\)
  1027. sed "s/^X//" >'OrderDither.c' <<'END_OF_FILE'
  1028. X/*
  1029. XOrdered Dithering
  1030. Xby Stephen Hawley
  1031. Xfrom "Graphics Gems", Academic Press, 1990
  1032. X*/
  1033. X
  1034. X/* Program to generate dithering matrices.
  1035. X * written by Jim Blandy, Oberlin College, jimb@occs.oberlin.edu
  1036. X * Gifted to, documented and revised by Stephen Hawley,
  1037. X * sdh@flash.bellcore.com
  1038. X * 
  1039. X * Generates a dithering matrix from the command line arguments.
  1040. X * The first argument, size, determines the dimensions of the 
  1041. X * matrix: 2^size by 2^size
  1042. X * The optional range argument is the range of values to be 
  1043. X * dithered over. By default, it is (2^size)^2, or simply the
  1044. X * total number of elements in the matrix.
  1045. X * The final output is suitable for inclusion in a C program.
  1046. X * A typical dithering function is something like this:
  1047. X * extern int dm[], size;
  1048. X *
  1049. X * int
  1050. X * dither(x,y, level)
  1051. X * register int x,y, level;
  1052. X * {
  1053. X *    return(level > dm[(x % size) + size * (y % size)];
  1054. X * }
  1055. X */
  1056. X
  1057. Xmain(argc, argv)
  1058. Xint argc;
  1059. Xchar **argv;
  1060. X{
  1061. X    register int size, range;
  1062. X
  1063. X    if (argc >= 2) size = atoi(argv[1]);
  1064. X    else size = 2;
  1065. X
  1066. X    if (argc == 3) range = atoi(argv[2]);
  1067. X    else range = (1 << size) * (1 << size);
  1068. X
  1069. X    printdither (size, range);
  1070. X}
  1071. X
  1072. X
  1073. Xprintdither (size, range)
  1074. Xregister int size, range;
  1075. X{
  1076. X    register int l = (1 << size), i;
  1077. X     /*
  1078. X      * print a dithering matrix.
  1079. X      * l is the length on a side.
  1080. X      */
  1081. X    range = range / (l * l);
  1082. X    puts("int dm[] = {");
  1083. X    for (i=0; i < l*l; i++) {
  1084. X        if (i % l == 0) /* tab in 4 spaces per row */
  1085. X            printf("    ");
  1086. X        /* print the dither value for this location
  1087. X          * scaled to the given range
  1088. X         */
  1089. X        printf("%4d", range * dithervalue(i / l, i % l, size));
  1090. X
  1091. X        /* commas after all but the last */
  1092. X        if (i + 1 < l * l)
  1093. X            putchar(',');
  1094. X        /* newline at the end of the row */
  1095. X        if ((i + 1) % l == 0)
  1096. X            putchar('\n');
  1097. X    }
  1098. X    puts("\n}; ");
  1099. X}
  1100. Xdithervalue(x, y, size)
  1101. Xregister int x, y, size;
  1102. X{
  1103. X    register int d;
  1104. X    /*
  1105. X     * calculate the dither value at a particular
  1106. X     * (x, y) over the size of the matrix.
  1107. X     */
  1108. X    d=0;
  1109. X    while (size-->0)    {
  1110. X        /* Think of d as the density. At every iteration,
  1111. X         * d is shifted left one and a new bit is put in the
  1112. X         * low bit based on x and y. If x is odd and y is even,
  1113. X         * or x is even and y is odd, a bit is put in. This
  1114. X         * generates the checkerboard seen in dithering.
  1115. X         * This quantity is shifted left again and the low bit of
  1116. X         * y is added in.
  1117. X         * This whole thing interleaves a checkerboard bit pattern
  1118. X         * and y's bits, which is the value you want.
  1119. X         */
  1120. X        d = (d <<1 | (x&1 ^ y&1))<<1 | y&1;
  1121. X        x >>= 1;
  1122. X        y >>= 1;
  1123. X    }
  1124. X    return(d);
  1125. X}
  1126. X
  1127. END_OF_FILE
  1128. if test 2483 -ne `wc -c <'OrderDither.c'`; then
  1129.     echo shar: \"'OrderDither.c'\" unpacked with wrong size!
  1130. fi
  1131. # end of 'OrderDither.c'
  1132. fi
  1133. if test -f 'PntOnLine.c' -a "${1}" != "-c" ; then 
  1134.   echo shar: Will not clobber existing file \"'PntOnLine.c'\"
  1135. else
  1136. echo shar: Extracting \"'PntOnLine.c'\" \(2024 characters\)
  1137. sed "s/^X//" >'PntOnLine.c' <<'END_OF_FILE'
  1138. X/* 
  1139. XA Fast 2D Point-On-Line Test
  1140. Xby Alan Paeth
  1141. Xfrom "Graphics Gems", Academic Press, 1990
  1142. X*/
  1143. X
  1144. X#include "GraphicsGems.h"
  1145. X
  1146. Xint PntOnLine(px,py,qx,qy,tx,ty)
  1147. X   int px, py, qx, qy, tx, ty;
  1148. X   {
  1149. X/*
  1150. X * given a line through P:(px,py) Q:(qx,qy) and T:(tx,ty)
  1151. X * return 0 if T is not on the line through      <--P--Q-->
  1152. X *        1 if T is on the open ray ending at P: <--P
  1153. X *        2 if T is on the closed interior along:   P--Q
  1154. X *        3 if T is on the open ray beginning at Q:    Q-->
  1155. X *
  1156. X * Example: consider the line P = (3,2), Q = (17,7). A plot
  1157. X * of the test points T(x,y) (with 0 mapped onto '.') yields:
  1158. X *
  1159. X *     8| . . . . . . . . . . . . . . . . . 3 3
  1160. X *  Y  7| . . . . . . . . . . . . . . 2 2 Q 3 3    Q = 2
  1161. X *     6| . . . . . . . . . . . 2 2 2 2 2 . . .
  1162. X *  a  5| . . . . . . . . 2 2 2 2 2 2 . . . . .
  1163. X *  x  4| . . . . . 2 2 2 2 2 2 . . . . . . . .
  1164. X *  i  3| . . . 2 2 2 2 2 . . . . . . . . . . .
  1165. X *  s  2| 1 1 P 2 2 . . . . . . . . . . . . . .    P = 2
  1166. X *     1| 1 1 . . . . . . . . . . . . . . . . .
  1167. X *      +--------------------------------------
  1168. X *        1 2 3 4 5 X-axis 10        15      19
  1169. X *
  1170. X * Point-Line distance is normalized with the Infinity Norm
  1171. X * avoiding square-root code and tightening the test vs the
  1172. X * Manhattan Norm. All math is done on the field of integers.
  1173. X * The latter replaces the initial ">= MAX(...)" test with
  1174. X * "> (ABS(qx-px) + ABS(qy-py))" loosening both inequality
  1175. X * and norm, yielding a broader target line for selection.
  1176. X * The tightest test is employed here for best discrimination
  1177. X * in merging collinear (to grid coordinates) vertex chains
  1178. X * into a larger, spanning vectors within the Lemming editor.
  1179. X */
  1180. X
  1181. X
  1182. X    if ( ABS((qy-py)*(tx-px)-(ty-py)*(qx-px)) >=
  1183. X        (MAX(ABS(qx-px), ABS(qy-py)))) return(0);
  1184. X    if (((qx<px)&&(px<tx)) || ((qy<py)&&(py<ty))) return(1);
  1185. X    if (((tx<px)&&(px<qx)) || ((ty<py)&&(py<qy))) return(1);
  1186. X    if (((px<qx)&&(qx<tx)) || ((py<qy)&&(qy<ty))) return(3);
  1187. X    if (((tx<qx)&&(qx<px)) || ((ty<qy)&&(qy<py))) return(3);
  1188. X    return(2);
  1189. X    }
  1190. END_OF_FILE
  1191. if test 2024 -ne `wc -c <'PntOnLine.c'`; then
  1192.     echo shar: \"'PntOnLine.c'\" unpacked with wrong size!
  1193. fi
  1194. # end of 'PntOnLine.c'
  1195. fi
  1196. if test -f 'PolyScan/fancytest.c' -a "${1}" != "-c" ; then 
  1197.   echo shar: Will not clobber existing file \"'PolyScan/fancytest.c'\"
  1198. else
  1199. echo shar: Extracting \"'PolyScan/fancytest.c'\" \(3298 characters\)
  1200. sed "s/^X//" >'PolyScan/fancytest.c' <<'END_OF_FILE'
  1201. X/*
  1202. X * fancytest.c: subroutine illustrating the use of poly_clip and poly_scan
  1203. X * for Phong-shading and texture mapping.
  1204. X *
  1205. X * Note: lines enclosed in angle brackets '<', '>' should be replaced
  1206. X * with the code described.
  1207. X * Makes calls to hypothetical packages "shade", "image", "texture", "zbuffer".
  1208. X *
  1209. X * Paul Heckbert    Dec 1989
  1210. X */
  1211. X
  1212. X#include <stdio.h>
  1213. X#include <math.h>
  1214. X#include "poly.h"
  1215. X
  1216. X#define XMAX 1280    /* hypothetical image width */
  1217. X#define YMAX 1024    /* hypothetical image height */
  1218. X#define LIGHT_INTEN 255        /* light source intensity */
  1219. X
  1220. Xvoid pixelproc();
  1221. X
  1222. Xfancytest()
  1223. X{
  1224. X    int i;
  1225. X    double WS[4][4];    /* world space to screen space transform */
  1226. X    Poly p;        /* a polygon */
  1227. X    Poly_vert *v;
  1228. X
  1229. X    static Poly_box box = {0, XMAX, 0, YMAX, -32678, 32767.99};
  1230. X    /* 3-D screen clipping box, continuous coordinates */
  1231. X
  1232. X    static Window win = {0, 0, XMAX-1, YMAX-1};
  1233. X    /* 2-D screen clipping window, discrete coordinates */
  1234. X
  1235. X    <initialize world space position (x,y,z), normal (nx,ny,nz), and texture
  1236. X     position (u,v) at each vertex of p; set p.n>
  1237. X    <set WS to world-to-screen transform>
  1238. X
  1239. X    /* transform vertices from world space to homogeneous screen space */
  1240. X    for (i=0; i<p.n; i++) {
  1241. X    v = &p.vert[i];
  1242. X    mx4_transform(v->x, v->y, v->z, 1., WS, &v->sx, &v->sy, &v->sz, &v->sw);
  1243. X    }
  1244. X
  1245. X    /* interpolate sx, sy, sz, sw, nx, ny, nz, u, v in poly_clip */
  1246. X    p.mask = POLY_MASK(sx) | POLY_MASK(sy) | POLY_MASK(sz) | POLY_MASK(sw) |
  1247. X    POLY_MASK(nx) | POLY_MASK(ny) | POLY_MASK(nz) |
  1248. X    POLY_MASK(u) | POLY_MASK(v);
  1249. X
  1250. X    poly_print("before clipping", &p);
  1251. X    if (poly_clip_to_box(&p, &box) == POLY_CLIP_OUT)    /* clip polygon */
  1252. X    return;                        /* quit if off-screen */
  1253. X
  1254. X    /* do homogeneous division of screen position and texture position */
  1255. X    for (i=0; i<p.n; i++) {
  1256. X    v = &p.vert[i];
  1257. X    v->sx /= v->sw;
  1258. X    v->sy /= v->sw;
  1259. X    v->sz /= v->sw;
  1260. X    v->u /= v->sw;
  1261. X    v->v /= v->sw;
  1262. X    v->q = 1./v->sw;
  1263. X    }
  1264. X    /*
  1265. X     * previously we ignored q (assumed q=1), henceforth ignore sw (assume sw=1)
  1266. X     * Interpolate sx, sy, sz, nx, ny, nz, u, v, q in poly_scan
  1267. X     */
  1268. X    p.mask &= ~POLY_MASK(sw);
  1269. X    p.mask |= POLY_MASK(q);
  1270. X
  1271. X    poly_print("scan converting the polygon", &p);
  1272. X    poly_scan(&p, &win, pixelproc);    /* scan convert! */
  1273. X}
  1274. X
  1275. Xstatic void pixelproc(x, y, pt)        /* called at each pixel by poly_scan */
  1276. Xint x, y;
  1277. XPoly_vert *pt;
  1278. X{
  1279. X    int sz, u, v, inten;
  1280. X    double len, nor[3], diff, spec;
  1281. X
  1282. X    sz = pt->sz;
  1283. X    if (sz < zbuffer_read(x, y)) {
  1284. X    len = sqrt(pt->nx*pt->nx + pt->ny*pt->ny + pt->nz*pt->nz);
  1285. X    nor[0] = pt->nx/len;        /* unitize the normal vector */
  1286. X    nor[1] = pt->ny/len;
  1287. X    nor[2] = pt->nz/len;
  1288. X    shade(nor, &diff, &spec);    /* compute specular and diffuse coeffs*/
  1289. X    u = pt->u/pt->q;        /* do homogeneous div. of texture pos */
  1290. X    v = pt->v/pt->q;
  1291. X    inten = texture_read(u, v)*diff + LIGHT_INTEN*spec;
  1292. X    image_write(x, y, inten);
  1293. X    zbuffer_write(x, y, sz);
  1294. X    }
  1295. X}
  1296. X
  1297. X/* mx4_transform: transform 4-vector p by matrix m yielding q: q = p*m */
  1298. X
  1299. Xmx4_transform(px, py, pz, pw, m, qxp, qyp, qzp, qwp)
  1300. Xdouble px, py, pz, pw, m[4][4], *qxp, *qyp, *qzp, *qwp;
  1301. X{
  1302. X    *qxp = px*m[0][0] + py*m[1][0] + pz*m[2][0] + pw*m[3][0];
  1303. X    *qyp = px*m[0][1] + py*m[1][1] + pz*m[2][1] + pw*m[3][1];
  1304. X    *qzp = px*m[0][2] + py*m[1][2] + pz*m[2][2] + pw*m[3][2];
  1305. X    *qwp = px*m[0][3] + py*m[1][3] + pz*m[2][3] + pw*m[3][3];
  1306. X}
  1307. END_OF_FILE
  1308. if test 3298 -ne `wc -c <'PolyScan/fancytest.c'`; then
  1309.     echo shar: \"'PolyScan/fancytest.c'\" unpacked with wrong size!
  1310. fi
  1311. # end of 'PolyScan/fancytest.c'
  1312. fi
  1313. if test -f 'PolyScan/poly.c' -a "${1}" != "-c" ; then 
  1314.   echo shar: Will not clobber existing file \"'PolyScan/poly.c'\"
  1315. else
  1316. echo shar: Extracting \"'PolyScan/poly.c'\" \(2310 characters\)
  1317. sed "s/^X//" >'PolyScan/poly.c' <<'END_OF_FILE'
  1318. X/*
  1319. X * poly.c: simple utilities for polygon data structures
  1320. X */
  1321. X
  1322. X#include "poly.h"
  1323. X
  1324. XPoly_vert *poly_dummy;        /* used superficially by POLY_MASK macro */
  1325. X
  1326. X/*
  1327. X * poly_print: print Poly p to stdout, prefixed by the label str
  1328. X */
  1329. X
  1330. Xvoid poly_print(str, p)
  1331. Xchar *str;
  1332. XPoly *p;
  1333. X{
  1334. X    int i;
  1335. X
  1336. X    printf("%s: %d sides\n", str, p->n);
  1337. X    poly_vert_label("        ", p->mask);
  1338. X    for (i=0; i<p->n; i++) {
  1339. X    printf("   v[%d] ", i);
  1340. X    poly_vert_print("", &p->vert[i], p->mask);
  1341. X    }
  1342. X}
  1343. X
  1344. Xvoid poly_vert_label(str, mask)
  1345. Xchar *str;
  1346. Xint mask;
  1347. X{
  1348. X    printf("%s", str);
  1349. X    if (mask&POLY_MASK(sx))   printf("   sx  ");
  1350. X    if (mask&POLY_MASK(sy))   printf("   sy  ");
  1351. X    if (mask&POLY_MASK(sz))   printf("   sz  ");
  1352. X    if (mask&POLY_MASK(sw))   printf("   sw  ");
  1353. X    if (mask&POLY_MASK(x))    printf("   x   ");
  1354. X    if (mask&POLY_MASK(y))    printf("   y   ");
  1355. X    if (mask&POLY_MASK(z))    printf("   z   ");
  1356. X    if (mask&POLY_MASK(u))    printf("   u   ");
  1357. X    if (mask&POLY_MASK(v))    printf("   v   ");
  1358. X    if (mask&POLY_MASK(q))    printf("   q   ");
  1359. X    if (mask&POLY_MASK(r))    printf("   r   ");
  1360. X    if (mask&POLY_MASK(g))    printf("   g   ");
  1361. X    if (mask&POLY_MASK(b))    printf("   b   ");
  1362. X    if (mask&POLY_MASK(nx))   printf("   nx  ");
  1363. X    if (mask&POLY_MASK(ny))   printf("   ny  ");
  1364. X    if (mask&POLY_MASK(nz))   printf("   nz  ");
  1365. X    printf("\n");
  1366. X}
  1367. X
  1368. Xvoid poly_vert_print(str, v, mask)
  1369. Xchar *str;
  1370. XPoly_vert *v;
  1371. Xint mask;
  1372. X{
  1373. X    printf("%s", str);
  1374. X    if (mask&POLY_MASK(sx)) printf(" %6.1f", v->sx);
  1375. X    if (mask&POLY_MASK(sy)) printf(" %6.1f", v->sy);
  1376. X    if (mask&POLY_MASK(sz)) printf(" %6.2f", v->sz);
  1377. X    if (mask&POLY_MASK(sw)) printf(" %6.2f", v->sw);
  1378. X    if (mask&POLY_MASK(x))  printf(" %6.2f", v->x);
  1379. X    if (mask&POLY_MASK(y))  printf(" %6.2f", v->y);
  1380. X    if (mask&POLY_MASK(z))  printf(" %6.2f", v->z);
  1381. X    if (mask&POLY_MASK(u))  printf(" %6.2f", v->u);
  1382. X    if (mask&POLY_MASK(v))  printf(" %6.2f", v->v);
  1383. X    if (mask&POLY_MASK(q))  printf(" %6.2f", v->q);
  1384. X    if (mask&POLY_MASK(r))  printf(" %6.4f", v->r);
  1385. X    if (mask&POLY_MASK(g))  printf(" %6.4f", v->g);
  1386. X    if (mask&POLY_MASK(b))  printf(" %6.4f", v->b);
  1387. X    if (mask&POLY_MASK(nx)) printf(" %6.3f", v->nx);
  1388. X    if (mask&POLY_MASK(ny)) printf(" %6.3f", v->ny);
  1389. X    if (mask&POLY_MASK(nz)) printf(" %6.3f", v->nz);
  1390. X    printf("\n");
  1391. X}
  1392. END_OF_FILE
  1393. if test 2310 -ne `wc -c <'PolyScan/poly.c'`; then
  1394.     echo shar: \"'PolyScan/poly.c'\" unpacked with wrong size!
  1395. fi
  1396. # end of 'PolyScan/poly.c'
  1397. fi
  1398. if test -f 'PolyScan/poly.h' -a "${1}" != "-c" ; then 
  1399.   echo shar: Will not clobber existing file \"'PolyScan/poly.h'\"
  1400. else
  1401. echo shar: Extracting \"'PolyScan/poly.h'\" \(1914 characters\)
  1402. sed "s/^X//" >'PolyScan/poly.h' <<'END_OF_FILE'
  1403. X/* poly.h: definitions for polygon package */
  1404. X
  1405. X#ifndef POLY_HDR
  1406. X#define POLY_HDR
  1407. X
  1408. X#define POLY_NMAX 8        /* max #sides to a polygon; change if needed */
  1409. X
  1410. Xtypedef struct {        /* A POLYGON VERTEX */
  1411. X    double sx, sy, sz, sw;    /* screen space position (sometimes homo.) */
  1412. X    double x, y, z;        /* world space position */
  1413. X    double u, v, q;        /* texture position (sometimes homogeneous) */
  1414. X    double r, g, b;        /* (red,green,blue) color */
  1415. X    double nx, ny, nz;        /* world space normal vector */
  1416. X} Poly_vert;
  1417. X/* update poly.c if you change this structure */
  1418. X
  1419. Xtypedef struct {        /* A POLYGON */
  1420. X    int n;            /* number of sides */
  1421. X    int mask;            /* interpolation mask for vertex elems */
  1422. X    Poly_vert vert[POLY_NMAX];    /* vertices */
  1423. X} Poly;
  1424. X/*
  1425. X * mask is an interpolation mask whose kth bit indicates whether the kth
  1426. X * double in a Poly_vert is relevant.
  1427. X * For example, if the valid attributes are sx, sy, and sz, then set
  1428. X *    mask = POLY_MASK(sx) | POLY_MASK(sy) | POLY_MASK(sz);
  1429. X */
  1430. X
  1431. Xtypedef struct {        /* A BOX (TYPICALLY IN SCREEN SPACE) */
  1432. X    double x0, x1;        /* left and right */
  1433. X    double y0, y1;        /* top and bottom */
  1434. X    double z0, z1;        /* near and far */
  1435. X} Poly_box;
  1436. X
  1437. Xtypedef struct {        /* WINDOW: A DISCRETE 2-D RECTANGLE */
  1438. X    int x0, y0;            /* xmin and ymin */
  1439. X    int x1, y1;            /* xmax and ymax (inclusive) */
  1440. X} Window;
  1441. X
  1442. X#define POLY_MASK(elem) (1 << (&poly_dummy->elem - (double *)poly_dummy))
  1443. X
  1444. X#define POLY_CLIP_OUT 0        /* polygon entirely outside box */
  1445. X#define POLY_CLIP_PARTIAL 1    /* polygon partially inside */
  1446. X#define POLY_CLIP_IN 2        /* polygon entirely inside box */
  1447. X
  1448. Xextern Poly_vert *poly_dummy;    /* used superficially by POLY_MASK macro */
  1449. X
  1450. Xvoid    poly_print(/* str, p */);
  1451. Xvoid    poly_vert_label(/* str, mask */);
  1452. Xvoid    poly_vert_print(/* str, v, mask */);
  1453. Xint    poly_clip_to_box(/* p1, box */);
  1454. Xvoid    poly_clip_to_halfspace(/* p, q, index, sign, k, name */);
  1455. Xvoid    poly_scan(/* p, win, pixelproc */);
  1456. X
  1457. X#endif
  1458. END_OF_FILE
  1459. if test 1914 -ne `wc -c <'PolyScan/poly.h'`; then
  1460.     echo shar: \"'PolyScan/poly.h'\" unpacked with wrong size!
  1461. fi
  1462. # end of 'PolyScan/poly.h'
  1463. fi
  1464. if test -f 'PolyScan/scantest.c' -a "${1}" != "-c" ; then 
  1465.   echo shar: Will not clobber existing file \"'PolyScan/scantest.c'\"
  1466. else
  1467. echo shar: Extracting \"'PolyScan/scantest.c'\" \(1854 characters\)
  1468. sed "s/^X//" >'PolyScan/scantest.c' <<'END_OF_FILE'
  1469. X/*
  1470. X * scantest.c: use poly_scan() for Gouraud shading and z-buffer demo.
  1471. X * Given the screen space X, Y, and Z of N-gon on command line,
  1472. X * print out all pixels during scan conversion.
  1473. X * This code could easily be modified to actually read and write pixels.
  1474. X *
  1475. X * Paul Heckbert    Dec 1989
  1476. X */
  1477. X
  1478. X#include <stdio.h>
  1479. X#include <math.h>
  1480. X#include "poly.h"
  1481. X
  1482. X#define XMAX 1280    /* hypothetical image width */
  1483. X#define YMAX 1024    /* hypothetical image height */
  1484. X
  1485. X#define FRANDOM() ((rand()&32767)/32767.)   /* random number between 0 and 1 */
  1486. X
  1487. Xstatic void pixelproc();
  1488. X
  1489. Xmain(ac, av)
  1490. Xint ac;
  1491. Xchar **av;
  1492. X{
  1493. X    int i;
  1494. X    Poly p;
  1495. X    static Window win = {0, 0, XMAX-1, YMAX-1};    /* screen clipping window */
  1496. X
  1497. X    if (ac<2 || ac != 2+3*(p.n = atoi(av[1]))) {
  1498. X    fprintf(stderr, "Usage: scantest N X1 Y1 Z1 X2 Y2 Z2 ... XN YN ZN\n");
  1499. X    exit(1);
  1500. X    }
  1501. X    for (i=0; i<p.n; i++) {
  1502. X    p.vert[i].sx = atof(av[2+3*i]);    /* set screen space x,y,z */
  1503. X    p.vert[i].sy = atof(av[3+3*i]);
  1504. X    p.vert[i].sz = atof(av[4+3*i]);
  1505. X    p.vert[i].r = FRANDOM();    /* random vertex colors, for kicks */
  1506. X    p.vert[i].g = FRANDOM();
  1507. X    p.vert[i].b = FRANDOM();
  1508. X    }
  1509. X    /* interpolate sx, sy, sz, r, g, and b in poly_scan */
  1510. X    p.mask = POLY_MASK(sx) | POLY_MASK(sy) | POLY_MASK(sz) |
  1511. X    POLY_MASK(r) | POLY_MASK(g) | POLY_MASK(b);
  1512. X
  1513. X    poly_print("scan converting the polygon", &p);
  1514. X    
  1515. X    poly_scan(&p, &win, pixelproc);    /* scan convert! */
  1516. X}
  1517. X
  1518. Xstatic void pixelproc(x, y, point)    /* called at each pixel by poly_scan */
  1519. Xint x, y;
  1520. XPoly_vert *point;
  1521. X{
  1522. X    printf("pixel (%d,%d) screenz=%g rgb=(%g,%g,%g)\n",
  1523. X    x, y, point->sz, point->r, point->g, point->b);
  1524. X
  1525. X    /*
  1526. X     * in real graphics program you could read and write pixels, e.g.:
  1527. X     *
  1528. X     *    if (point->sz < zbuffer_read(x, y)) {
  1529. X     *        image_write_rgb(x, y, point->r, point->g, point->b);
  1530. X     *        zbuffer_write(x, y, point->sz);
  1531. X     *  }
  1532. X     */
  1533. X}
  1534. END_OF_FILE
  1535. if test 1854 -ne `wc -c <'PolyScan/scantest.c'`; then
  1536.     echo shar: \"'PolyScan/scantest.c'\" unpacked with wrong size!
  1537. fi
  1538. # end of 'PolyScan/scantest.c'
  1539. fi
  1540. if test -f 'Quaternions.c' -a "${1}" != "-c" ; then 
  1541.   echo shar: Will not clobber existing file \"'Quaternions.c'\"
  1542. else
  1543. echo shar: Extracting \"'Quaternions.c'\" \(2216 characters\)
  1544. sed "s/^X//" >'Quaternions.c' <<'END_OF_FILE'
  1545. X/*
  1546. XUsing Quaternions for Coding 3D Transformations
  1547. XPatrick-Gilles Maillot
  1548. Xfrom "Graphics Gems", Academic Press, 1990
  1549. X*/
  1550. X
  1551. Xextern double P[3], Q[4], M[4][4];
  1552. X
  1553. Xset_obs_position(x,y,z)
  1554. Xfloat    x, y, z;
  1555. X{
  1556. Xint    i;
  1557. X
  1558. X/*
  1559. X * Set the values of the eye's position.
  1560. X * The position here represents the position of the orthonormal base
  1561. X * in respect to the observer.
  1562. X */
  1563. X       P[0] = -x;
  1564. X       P[1] = -y;
  1565. X       P[2] = -z;
  1566. X/*
  1567. X * Set the visualization to be in the decreasing x axis
  1568. X */
  1569. X    Q[0] = 1.;
  1570. X      for (i = 1; i < 4; i++) Q[i] = 0.;
  1571. X}
  1572. X
  1573. Xtranslate_quaternion(x,i,w)
  1574. Xfloat    x;
  1575. Xint    i, w;
  1576. X{
  1577. Xint    j, k;
  1578. Xfloat    A, B, D, E, F;
  1579. X    
  1580. X    if (w < 0) {
  1581. X/*
  1582. X * The observer moves in respect to the scene.
  1583. X */
  1584. X    P[i - 1] -= x;
  1585. X  } else {
  1586. X
  1587. X/* 
  1588. X * The scene moves in respect to the observer.
  1589. X * Compute the successor axis of i [1,2,3];
  1590. X * and then the successor axis of j [1,2,3];
  1591. X */
  1592. X    if ((j = i + 1) > 3) j = 1;
  1593. X    if ((k = j + 1) > 3) k = 1;
  1594. X    A = Q[j]; B = Q[k]; F = Q[0]; E = Q[i];
  1595. X    P[i - 1] += x * (E * E + F * F - A * A - B * B);
  1596. X    D = x + x;
  1597. X    P[j - 1] += D * (E * A + F * B);
  1598. X    P[k - 1] += D * (E * B + F * A);
  1599. X  }
  1600. X}
  1601. X
  1602. Xrotate_quaternion(x,y,i,w)
  1603. Xfloat    x, y;
  1604. Xint    i, w;
  1605. X{
  1606. Xint    j, k;
  1607. Xfloat    E, F, R1;
  1608. X/*
  1609. X * Compute the successor axis of i [1,2,3] and  j [1,2,3];
  1610. X */
  1611. X    if ((j = i + 1) > 3) j = 1;
  1612. X    if ((k = j + 1) > 3) k = 1;
  1613. X    E = Q[i];
  1614. X    Q[i] = E * x + w * y * Q[0];
  1615. X    Q[0] = Q[0] * x - w * y * E;
  1616. X    E = Q[j];
  1617. X    Q[j] = E * x + y * Q[k];
  1618. X    Q[k] = Q[k] * x - y * E;
  1619. X      if (w < 0) {
  1620. X/* Compute a new position if the observer moves in respect to the scene. */
  1621. X        j -= 1; k -= 1;
  1622. X        R1 = x * x - y * y;
  1623. X        F = 2. * x * y;
  1624. X        E = P[j];
  1625. X        P[j] = E * R1 + F * P[k];
  1626. X        P[k] = P[k] * R1 - F * E;
  1627. X      }
  1628. X}
  1629. X
  1630. X
  1631. XEvaluate_matrix()
  1632. X{
  1633. Xfloat    e, f, r[4];
  1634. Xint    i, j, k;
  1635. X/*
  1636. X * We will need some square values!
  1637. X */
  1638. X    for (i = 0; i < 4; i++) r[i] = Q[i] * Q[i];
  1639. X/*
  1640. X * Compute each element of the matrix.
  1641. X * j is the successor of i (in 1,2,3), while k is the successor of j.
  1642. X */
  1643. X      for (i = 1; i < 4; i++) {
  1644. X        if ((j = i + 1) > 3) j = 1;
  1645. X        if ((k = j + 1) > 3) k = 1;
  1646. X        e = 2. * Q[i] * Q[j];
  1647. X        f = 2. * Q[k] * Q[0];
  1648. X        M[j][i] = e - f;
  1649. X        M[i][j] = e + f;
  1650. X        M[i][i] = r[i] + r[0] - r[j] - r[k];
  1651. X        M[0][i] = P[i - 1];
  1652. X        M[i][0] = 0.;
  1653. X      }
  1654. X    M[0][0] = 1.;
  1655. X}
  1656. X
  1657. X
  1658. X
  1659. X
  1660. END_OF_FILE
  1661. if test 2216 -ne `wc -c <'Quaternions.c'`; then
  1662.     echo shar: \"'Quaternions.c'\" unpacked with wrong size!
  1663. fi
  1664. # end of 'Quaternions.c'
  1665. fi
  1666. if test -f 'SeedFill.c' -a "${1}" != "-c" ; then 
  1667.   echo shar: Will not clobber existing file \"'SeedFill.c'\"
  1668. else
  1669. echo shar: Extracting \"'SeedFill.c'\" \(2371 characters\)
  1670. sed "s/^X//" >'SeedFill.c' <<'END_OF_FILE'
  1671. X/*
  1672. X * A Seed Fill Algorithm
  1673. X * by Paul Heckbert
  1674. X * from "Graphics Gems", Academic Press, 1990
  1675. X */
  1676. X
  1677. X/*
  1678. X * fill.c : simple seed fill program
  1679. X * Calls pixelread() to read pixels, pixelwrite() to write pixels.
  1680. X *
  1681. X * Paul Heckbert    13 Sept 1982, 28 Jan 1987
  1682. X */
  1683. X
  1684. Xtypedef struct {        /* window: a discrete 2-D rectangle */
  1685. X    int x0, y0;            /* xmin and ymin */
  1686. X    int x1, y1;            /* xmax and ymax (inclusive) */
  1687. X} Window;
  1688. X
  1689. Xtypedef int Pixel;        /* 1-channel frame buffer assumed */
  1690. X
  1691. XPixel pixelread();
  1692. X
  1693. Xtypedef struct {short y, xl, xr, dy;} Segment;
  1694. X/*
  1695. X * Filled horizontal segment of scanline y for xl<=x<=xr.
  1696. X * Parent segment was on line y-dy.  dy=1 or -1
  1697. X */
  1698. X
  1699. X#define MAX 10000        /* max depth of stack */
  1700. X
  1701. X#define PUSH(Y, XL, XR, DY)    /* push new segment on stack */ \
  1702. X    if (sp<stack+MAX && Y+(DY)>=win->y0 && Y+(DY)<=win->y1) \
  1703. X    {sp->y = Y; sp->xl = XL; sp->xr = XR; sp->dy = DY; sp++;}
  1704. X
  1705. X#define POP(Y, XL, XR, DY)    /* pop segment off stack */ \
  1706. X    {sp--; Y = sp->y+(DY = sp->dy); XL = sp->xl; XR = sp->xr;}
  1707. X
  1708. X/*
  1709. X * fill: set the pixel at (x,y) and all of its 4-connected neighbors
  1710. X * with the same pixel value to the new pixel value nv.
  1711. X * A 4-connected neighbor is a pixel above, below, left, or right of a pixel.
  1712. X */
  1713. X
  1714. Xfill(x, y, win, nv)
  1715. Xint x, y;    /* seed point */
  1716. XWindow *win;    /* screen window */
  1717. XPixel nv;    /* new pixel value */
  1718. X{
  1719. X    int l, x1, x2, dy;
  1720. X    Pixel ov;    /* old pixel value */
  1721. X    Segment stack[MAX], *sp = stack;    /* stack of filled segments */
  1722. X
  1723. X    ov = pixelread(x, y);        /* read pv at seed point */
  1724. X    if (ov==nv || x<win->x0 || x>win->x1 || y<win->y0 || y>win->y1) return;
  1725. X    PUSH(y, x, x, 1);            /* needed in some cases */
  1726. X    PUSH(y+1, x, x, -1);        /* seed segment (popped 1st) */
  1727. X
  1728. X    while (sp>stack) {
  1729. X    /* pop segment off stack and fill a neighboring scan line */
  1730. X    POP(y, x1, x2, dy);
  1731. X    /*
  1732. X     * segment of scan line y-dy for x1<=x<=x2 was previously filled,
  1733. X     * now explore adjacent pixels in scan line y
  1734. X     */
  1735. X    for (x=x1; x>=win->x0 && pixelread(x, y)==ov; x--)
  1736. X        pixelwrite(x, y, nv);
  1737. X    if (x>=x1) goto skip;
  1738. X    l = x+1;
  1739. X    if (l<x1) PUSH(y, l, x1-1, -dy);        /* leak on left? */
  1740. X    x = x1+1;
  1741. X    do {
  1742. X        for (; x<=win->x1 && pixelread(x, y)==ov; x++)
  1743. X        pixelwrite(x, y, nv);
  1744. X        PUSH(y, l, x-1, dy);
  1745. X        if (x>x2+1) PUSH(y, x2+1, x-1, -dy);    /* leak on right? */
  1746. Xskip:        for (x++; x<=x2 && pixelread(x, y)!=ov; x++);
  1747. X        l = x;
  1748. X    } while (x<=x2);
  1749. X    }
  1750. X}
  1751. END_OF_FILE
  1752. if test 2371 -ne `wc -c <'SeedFill.c'`; then
  1753.     echo shar: \"'SeedFill.c'\" unpacked with wrong size!
  1754. fi
  1755. # end of 'SeedFill.c'
  1756. fi
  1757. if test -f 'SquareRoot.c' -a "${1}" != "-c" ; then 
  1758.   echo shar: Will not clobber existing file \"'SquareRoot.c'\"
  1759. else
  1760. echo shar: Extracting \"'SquareRoot.c'\" \(2048 characters\)
  1761. sed "s/^X//" >'SquareRoot.c' <<'END_OF_FILE'
  1762. X/*
  1763. X * A High Speed, Low Precision Square Root
  1764. X * by Paul Lalonde and Robert Dawson
  1765. X * from "Graphics Gems", Academic Press, 1990
  1766. X */
  1767. X
  1768. X
  1769. X/*
  1770. X * SPARC implementation of a fast square root by table 
  1771. X * lookup.
  1772. X * SPARC floating point format is as follows:
  1773. X *
  1774. X * BIT 31     30     23     22     0
  1775. X *     sign    exponent    mantissa
  1776. X */
  1777. X
  1778. X#include <math.h>
  1779. X
  1780. Xstatic short sqrttab[0x100];    /* declare table of square roots */
  1781. X
  1782. Xvoid
  1783. Xbuild_table()
  1784. X{
  1785. X    unsigned short i;
  1786. X    float f;
  1787. X    unsigned int *fi = &f;    /* to access the bits of a float in  */
  1788. X                /* C quickly we must misuse pointers */
  1789. X
  1790. X    for(i=0; i<= 0x7f; i++) {
  1791. X        *fi = 0;
  1792. X
  1793. X        /*
  1794. X         * Build a float with the bit pattern i as mantissa
  1795. X         * and an exponent of 0, stored as 127
  1796. X           */
  1797. X
  1798. X        *fi = (i << 16) | (127 << 23);
  1799. X        f = sqrt(f);
  1800. X
  1801. X        /*
  1802. X         * Take the square root then strip the first 7 bits of
  1803. X         * the mantissa into the table
  1804. X         */
  1805. X
  1806. X        sqrttab[i] = (*fi & 0x7fffff) >> 16;
  1807. X
  1808. X        /*
  1809. X         * Repeat the process, this time with an exponent of
  1810. X         * 1, stored as 128
  1811. X         */
  1812. X
  1813. X        *fi = 0;
  1814. X        *fi = (i << 16) | (128 << 23);
  1815. X        f = sqrt(f);
  1816. X        sqrttab[i+0x80] = (*fi & 0x7fffff) >> 16;
  1817. X    }
  1818. X}
  1819. X
  1820. X/*
  1821. X * fsqrt - fast square root by table lookup
  1822. X */
  1823. Xfloat
  1824. Xfsqrt(float n)
  1825. X{
  1826. X    unsigned int *num = &n;    /* to access the bits of a float in C
  1827. X                 * we must misuse pointers */
  1828. X                            
  1829. X    short e;        /* the exponent */
  1830. X    if (n == 0) return (0);    /* check for square root of 0 */
  1831. X    e = (*num >> 23) - 127;    /* get the exponent - on a SPARC the */
  1832. X                /* exponent is stored with 127 added */
  1833. X    *num & = 0x7fffff;    /* leave only the mantissa */
  1834. X    if (e & 0x01) *num | = 0x800000;
  1835. X                /* the exponent is odd so we have to */
  1836. X                /* look it up in the second half of  */
  1837. X                /* the lookup table, so we set the   */
  1838. X                /* high bit                   */
  1839. X    e >>= 1;        /* divide the exponent by two */
  1840. X                /* note that in C the shift */
  1841. X                /* operators are sign preserving */
  1842. X                /* for signed operands */
  1843. X    /* Do the table lookup, based on the quaternary mantissa,
  1844. X     * then reconstruct the result back into a float
  1845. X     */
  1846. X    *num = ((sqrttab[*num >> 16]) << 16) | ((e + 127) << 23);
  1847. X    return(n);
  1848. X}
  1849. END_OF_FILE
  1850. if test 2048 -ne `wc -c <'SquareRoot.c'`; then
  1851.     echo shar: \"'SquareRoot.c'\" unpacked with wrong size!
  1852. fi
  1853. # end of 'SquareRoot.c'
  1854. fi
  1855. if test -f 'Sturm/main.c' -a "${1}" != "-c" ; then 
  1856.   echo shar: Will not clobber existing file \"'Sturm/main.c'\"
  1857. else
  1858. echo shar: Extracting \"'Sturm/main.c'\" \(2173 characters\)
  1859. sed "s/^X//" >'Sturm/main.c' <<'END_OF_FILE'
  1860. X/*
  1861. XUsing Sturm Sequences to Bracket Real Roots of Polynomial Equations
  1862. Xby D.G. Hook and P.R. McAree
  1863. Xfrom "Graphics Gems", Academic Press, 1990
  1864. X*/
  1865. X
  1866. X/*
  1867. X * main.c
  1868. X *
  1869. X *    a sample driver program.
  1870. X */
  1871. X#include <stdio.h>
  1872. X#include <math.h>
  1873. X#include "solve.h"
  1874. X
  1875. X/*
  1876. X * a driver program for a root solver.
  1877. X */
  1878. Xmain()
  1879. X{
  1880. X    poly    sseq[MAX_ORDER];
  1881. X    double     min, max, roots[MAX_ORDER];
  1882. X    int        i, j, order, nroots, nchanges, np, atmin, atmax;
  1883. X
  1884. X    /*
  1885. X     * get the details...
  1886. X     */
  1887. X
  1888. X    printf("Please enter order of polynomial: ");
  1889. X    scanf("%d", &order);
  1890. X
  1891. X    printf("\n");
  1892. X
  1893. X    for (i = order; i >= 0; i--) {
  1894. X            printf("Please enter coefficient number %d: ", i);
  1895. X            scanf("%lf", &sseq[0].coef[i]);
  1896. X    }
  1897. X
  1898. X    printf("\n");
  1899. X
  1900. X    /*
  1901. X     * build the Sturm sequence
  1902. X     */
  1903. X    np = buildsturm(order, sseq);
  1904. X
  1905. X    printf("Sturm sequence for:\n");
  1906. X
  1907. X    for (i = order; i >= 0; i--)
  1908. X            printf("%lf ", sseq[0].coef[i]);
  1909. X
  1910. X    printf("\n\n");
  1911. X
  1912. X    for (i = 0; i <= np; i++) {
  1913. X            for (j = sseq[i].ord; j >= 0; j--)
  1914. X                printf("%lf ", sseq[i].coef[j]);
  1915. X            printf("\n");
  1916. X    }
  1917. X
  1918. X    printf("\n");
  1919. X
  1920. X
  1921. X    /* 
  1922. X     * get the number of real roots
  1923. X     */
  1924. X    nroots = numroots(np, sseq, &atmin, &atmax);
  1925. X
  1926. X    if (nroots == 0) {
  1927. X            printf("solve: no real roots\n");
  1928. X            exit(0);
  1929. X    }
  1930. X
  1931. X    /*
  1932. X     * calculate the bracket that the roots live in
  1933. X     */
  1934. X    min = -1.0;
  1935. X    nchanges = numchanges(np, sseq, min);
  1936. X    for (i = 0; nchanges != atmin && i != MAXPOW; i++) { 
  1937. X            min *= 10.0;
  1938. X            nchanges = numchanges(np, sseq, min);
  1939. X    }
  1940. X
  1941. X    if (nchanges != atmin) {
  1942. X            printf("solve: unable to bracket all negetive roots\n");
  1943. X            atmin = nchanges;
  1944. X    }
  1945. X
  1946. X    max = 1.0;
  1947. X    nchanges = numchanges(np, sseq, max);
  1948. X    for (i = 0; nchanges != atmax && i != MAXPOW; i++) { 
  1949. X            max *= 10.0;
  1950. X            nchanges = numchanges(np, sseq, max);
  1951. X    }
  1952. X
  1953. X    if (nchanges != atmax) {
  1954. X            printf("solve: unable to bracket all positive roots\n");
  1955. X            atmax = nchanges;
  1956. X    }
  1957. X
  1958. X    nroots = atmin - atmax;
  1959. X
  1960. X    /*
  1961. X     * perform the bisection.
  1962. X     */
  1963. X    sbisect(np, sseq, min, max, atmin, atmax, roots);
  1964. X
  1965. X
  1966. X    /*
  1967. X     * write out the roots...
  1968. X     */
  1969. X    if (nroots == 1) {
  1970. X            printf("\n1 distinct real root at x = %f\n", roots[0]);
  1971. X    } else {
  1972. X            printf("\n%d distinct real roots for x: ", nroots);
  1973. X
  1974. X            for (i = 0; i != nroots; i++)
  1975. X                printf("%f ", roots[i]);
  1976. X            printf("\n");
  1977. X    }
  1978. X}
  1979. END_OF_FILE
  1980. if test 2173 -ne `wc -c <'Sturm/main.c'`; then
  1981.     echo shar: \"'Sturm/main.c'\" unpacked with wrong size!
  1982. fi
  1983. # end of 'Sturm/main.c'
  1984. fi
  1985. if test -f 'TriPoints.c' -a "${1}" != "-c" ; then 
  1986.   echo shar: Will not clobber existing file \"'TriPoints.c'\"
  1987. else
  1988. echo shar: Extracting \"'TriPoints.c'\" \(2746 characters\)
  1989. sed "s/^X//" >'TriPoints.c' <<'END_OF_FILE'
  1990. X/* 
  1991. XGenerating Random Points in Triangles
  1992. Xby Greg Turk
  1993. Xfrom "Graphics Gems", Academic Press, 1990
  1994. X*/
  1995. X
  1996. X#include <math.h>
  1997. X#include "GraphicsGems.h"
  1998. X
  1999. X/*****************************************************************
  2000. XCompute relative areas of sub-triangles that form a convex polygon.
  2001. XThere are vcount-2 sub-triangles, each defined by the first point
  2002. Xin the polygon and two other adjacent points.
  2003. X
  2004. XThis procedure should be called once before using 
  2005. Xsquare_to_polygon().
  2006. X
  2007. XEntry:
  2008. X  vertices - list of the vertices of a convex polygon
  2009. X  vcount   - number of vertices of polygon
  2010. XExit:
  2011. X  areas - relative areas of sub-triangles of polygon
  2012. X*******************************************************************/
  2013. X
  2014. Xtriangle_areas (vertices, vcount, areas)
  2015. X      Point3 vertices[];
  2016. X    int vcount;
  2017. X      float areas[];
  2018. X{
  2019. X      int i;
  2020. X      float area_sum = 0;
  2021. X      Vector3 v1,v2,v3;
  2022. X
  2023. X  /* compute relative areas of the sub-triangles of polygon */
  2024. X
  2025. X      for (i = 0; i < vcount - 2; i++) {
  2026. X         V3Sub(&vertices[i+1], &vertices[0], &v1);
  2027. X         V3Sub(&vertices[i+2], &vertices[0], &v2);
  2028. X         V3Cross(&v1, &v2, &v3);
  2029. X         areas[i] = LengthVector3(&v3);
  2030. X         area_sum += areas[i];
  2031. X  }
  2032. X
  2033. X  /* normalize areas so that the sum of all sub-triangles is one */
  2034. X
  2035. X      for (i = 0; i < vcount - 2; i++)
  2036. X        areas[i] /= area_sum;
  2037. X}
  2038. X
  2039. X
  2040. X/******************************************************************
  2041. XMap a point from the square [0,1] x [0,1] into a convex polygon.
  2042. XUniform random points in the square will generate uniform random 
  2043. Xpoints in the polygon.
  2044. X
  2045. XThe procedure triangle_areas() must be called once to compute 
  2046. X'areas', and then this procedure can be called repeatedly.
  2047. X
  2048. XEntry:
  2049. X  vertices - list of the vertices of a convex polygon
  2050. X  vcount   - number of vertices of polygon
  2051. X  areas    - relative areas of sub-triangles of polygon
  2052. X  s,t      - position in the square [0,1] x [0,1]
  2053. XExit:
  2054. X  p - position in polygon
  2055. X*********************************************************************/
  2056. X
  2057. Xsquare_to_polygon (vertices, vcount, areas, s, t, p)
  2058. X     Point3 vertices[];
  2059. X      int vcount;
  2060. X      float areas[];
  2061. X      float s,t;
  2062. X      Point3 *p;
  2063. X{ 
  2064. X    int i;
  2065. X    float area_sum = 0;
  2066. X      float a,b,c;
  2067. X
  2068. X  /* use 's' to pick one sub-triangle, weighted by relative */
  2069. X  /* area of triangles */
  2070. X
  2071. X      for (i = 0; i < vcount - 2; i++) {
  2072. X        area_sum += areas[i];
  2073. X        if (area_sum >= s)
  2074. X              break;
  2075. X  }
  2076. X
  2077. X  /* map 's' into the interval [0,1] */
  2078. X
  2079. X  s = (s - area_sum + areas[i]) / areas[i];
  2080. X
  2081. X  /* map (s,t) to a point in that sub-triangle */
  2082. X
  2083. X  t = sqrt(t);
  2084. X
  2085. X  a = 1 - t;
  2086. X  b = (1 - s) * t;
  2087. X  c = s * t;
  2088. X
  2089. X  p->x = a * vertices[0].x + b * vertices[i+1].x + c * vertices[i+2].x;
  2090. X  p->y = a * vertices[0].y + b * vertices[i+1].y + c * vertices[i+2].y;
  2091. X  p->z = a * vertices[0].z + b * vertices[i+1].z + c * vertices[i+2].z;
  2092. X}
  2093. END_OF_FILE
  2094. if test 2746 -ne `wc -c <'TriPoints.c'`; then
  2095.     echo shar: \"'TriPoints.c'\" unpacked with wrong size!
  2096. fi
  2097. # end of 'TriPoints.c'
  2098. fi
  2099. echo shar: End of archive 2 \(of 5\).
  2100. cp /dev/null ark2isdone
  2101. MISSING=""
  2102. for I in 1 2 3 4 5 ; do
  2103.     if test ! -f ark${I}isdone ; then
  2104.     MISSING="${MISSING} ${I}"
  2105.     fi
  2106. done
  2107. if test "${MISSING}" = "" ; then
  2108.     echo You have unpacked all 5 archives.
  2109.     rm -f ark[1-9]isdone
  2110. else
  2111.     echo You still need to unpack the following archives:
  2112.     echo "        " ${MISSING}
  2113. fi
  2114. ##  End of shell archive.
  2115. exit 0
  2116.  
  2117.